# -*- coding: utf-8 -*-
"""
The python script to generate Figure 5. 
It includes 1. read monthly air temperature, huv, sensible heat flux, and 
latent heat flux data on three pairs of land tiles
2. Calculate the Bowen ratio from sensible and latent heat flux
3. Calculate the monthly air temperature, huv, and Bowen ratio in the 
North hemisphere for three pairs of land tiles
3. Plot the seasonal variations of air temperature, huv, and Bowen ratio for
three pairs of land tiles 
"""
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

Dir ='G:\\SSP5UC_DataPreparation\\deposit\\Hourly\\' # The directory you are working with
LonNum=288 # Number of longitudinal grids
LatNum=192 # Number of latitudinal grids

#Read the area of each grid in the North hemisphere
AreaFile=nc.Dataset('G:\\SSP5UC_DataPreparation\\deposit\\Surface\\b.e21.SSP585UC_IndividualSoil_surfdata_0.9x1.25_16pfts.nc')
area=AreaFile.variables['AREA'][96:192,:]

def AreaWeightedMean(VarData):
    """
    AreaWeightedMean
    Description: a function to calculate the area-weighed mean value of a variable
    Input： VarData--a 2D data in the shape of (96,288)
    Output: GlobalMean--area-weighed global mean value of VarData
     
    """  
    MaskedArea=np.ma.masked_where(np.ma.getmask(VarData[0]),area)
    GlobalMean=np.nansum(VarData*(MaskedArea/np.sum(MaskedArea)),axis=(1,2))
    return GlobalMean

def CalMeanDuirnalCycle(Var):
    """
    CalMeanDuirnalCycle
    Description: a function to calculate the long term mean duirnal cycles
    Input： Var--continuous hourly data during a certain time period in the shape of (hours,LatNum,LonNum)
    Output: VarMean--long-term mean duirnal cycles in the shape of (24,LatNum,LonNum)
     
    """     
    DayCount=int(len(Var[:,0,0])/24) # How many days 
    VarMean=np.empty((24,LatNum,LonNum));VarMean[:] = np.nan; VarMean= np.ma.masked_invalid(VarMean)
    for i in range(0,24): # loop through 0:00-23:00 and calculate the long-term mean values
        VarMean[i,:,:]=(np.nanmean(Var[i:24*(DayCount-1)+i+1:24,:,:],axis=0))    
    return VarMean


def ConvertUTC2LTC(UTCData):
    """ 
    Description: a function to convert data using UTC as data using Local Time
    Input: UTCData--the data in the shape of (24,LatNum,LonNum) using UTC
    Output: LTData--the data in the shape of (24,LatNum,LonNum) using LT
    """ 
    # Here TZLon is a 24*12 helper array. The values(0~287) in TZLon stands for the index on the longitudinal dimension.
    # The grids at the same row indicate longitudinal indexes belonging to the same Timezone. 
    # Their Timezone=row number when row number<12, and Timezone=row numbers-24 when row number>12
    TZLon=np.empty((24, int(LonNum/24)))
    TZInterval=int(LonNum/24) # There are 24 time zones. Every time zone contains LatNum*12 grids
    HalfInterval=int(TZInterval/2)
    for j in range(0,24):
        TZLon[j,:]=np.arange(int(HalfInterval+1+TZInterval*(j-1)),int(HalfInterval+1+TZInterval*(j))) -1 
    TZLon[0,0:HalfInterval]=np.arange(LonNum+1-HalfInterval,LonNum+1)-1
    TZLon=np.int16(TZLon)  
    
    # Initialize LTData as a 24*192*288 array
    LTData=np.empty((24,LatNum,LonNum));LTData[:] = np.nan; LTData= np.ma.masked_invalid(LTData)
    
    # For grids in each time zone, reorder the data such that the data at 12 AM local time is the first, and them 1AM, 2AM...
    for k in range(0,24):
        if k==0: # At time zone zero, make no adjustments
            LTData[0:24,:,LonNum-HalfInterval:LonNum]=UTCData[0:24,:,LonNum-HalfInterval:LonNum]
            LTData[0:24,:,0:HalfInterval]=UTCData[0:24,:,0:HalfInterval]
        elif k in range(1,24): # At time zone zero, make no adjustments
            Temp=24-k
            LTData[0:24,:,TZLon[k,0]:(TZLon[k,TZInterval-1]+1)]=np.ma.concatenate((UTCData[Temp:24,:,TZLon[k,0]:(TZLon[k,TZInterval-1]+1)],UTCData[0:Temp,:,TZLon[k,0]:(TZLon[k,TZInterval-1]+1)]),axis=0)

    # Apply the lake and rural masks just to make sure only grids with both lake and rural tiles are used
    LTData=np.ma.masked_where(np.ma.getmask(MaskArray),LTData) 
    return LTData

def ReadVar(Varname,Perturbed,Base):
    """
    ReadVar
    Description: a function to read hourly state varibles for the perturbed and base land tiles
                 and then calculate the duirnal cycles of their mean values in North hemisphere(NH) 
    Input： Varname--the short for the targeted variable (tas/huv)
    Perturbed--the name of the perturned land tiles (u_Urban/l_Lake/g_Grass)
    Base--the name of the base land tiles (r_Rural/r_Rural/t_Tree)
    Output: PerturbedMean--duirnal cycles of the NH mean values of the targeted variable on the perturbed land tile
            BaseMean--duirnal cycles of the NH mean values of the targeted variable on the base land tile
    """ 
    for year in range(2019,2023):
        Perturbedfile= nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_'+Varname+Perturbed+'_Hourly_'+str(year)+'.nc')
        tasp=Perturbedfile.variables[Varname+Perturbed[0]][:,:,:]
        Basefile= nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_'+Varname+Base+'_Hourly_'+str(year)+'.nc')
        tasb=Basefile.variables[Varname+Base[0]][:,:,:]
        # Calculate long-term duirnal cycles of this variable
        taspMean=CalMeanDuirnalCycle(tasp);
        tasbMean=CalMeanDuirnalCycle(tasb);
        # Convert UTC to Local Time
        taspMeanCon=ConvertUTC2LTC(taspMean)
        tasbMeanCon=ConvertUTC2LTC(tasbMean)
        # Add an addtional dimension so that the data can be stacked on this dimension
        taspMeanCon=taspMeanCon[np.newaxis,:,:,:]
        tasbMeanCon=tasbMeanCon[np.newaxis,:,:,:]
        if year==2019:
            PerturbedAll=taspMeanCon
            BaseAll=tasbMeanCon
        else:
            PerturbedAll=np.ma.concatenate((PerturbedAll,taspMeanCon),axis=0)
            BaseAll=np.ma.concatenate((BaseAll,tasbMeanCon),axis=0)
    # Five year mean
    PerturbedMean0=np.nanmean(PerturbedAll,axis=0)
    BaseMean0=np.nanmean(BaseAll,axis=0)
    # Here 96:192 means only use data in NH
    PerturbedMean=PerturbedMean0[:,96:192,:]
    BaseMean=BaseMean0[:,96:192,:]
    PerturbedMean=AreaWeightedMean(PerturbedMean)
    BaseMean=AreaWeightedMean(BaseMean)
    return PerturbedMean,BaseMean

def ReadVarNoMean(Varname,Perturbed,Base):
    """
    ReadVar
    Description: a function to read hourly state varibles for the perturbed and base land tiles
                 and then calculate the duirnal cycles of the NH result of the targeted variable
    Input： Varname--the short for the targeted variable (tas/huv)
    Perturbed--the name of the perturned land tiles (u_Urban/l_Lake/g_Grass)
    Base--the name of the base land tiles (r_Rural/r_Rural/t_Tree)
    Output: PerturbedMean--duirnal cycles of the NH result of the targeted variable on the perturbed land tile
            BaseMean--duirnal cycles of the NH result of the targeted variable on the base land tile
    """ 
    for year in range(2019,2023):
        Perturbedfile= nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_'+Varname+Perturbed+'_Hourly_'+str(year)+'.nc')
        tasp=Perturbedfile.variables[Varname+Perturbed[0]][:,:,:]
        Basefile= nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_'+Varname+Base+'_Hourly_'+str(year)+'.nc')
        tasb=Basefile.variables[Varname+Base[0]][:,:,:]
        # Calculate long-term duirnal cycles of this variable
        taspMean=CalMeanDuirnalCycle(tasp);
        tasbMean=CalMeanDuirnalCycle(tasb);
        # Convert UTC to Local Time
        taspMeanCon=ConvertUTC2LTC(taspMean)
        tasbMeanCon=ConvertUTC2LTC(tasbMean)
        # Add an addtional dimension so that the data can be stacked on this dimension
        taspMeanCon=taspMeanCon[np.newaxis,:,:,:]
        tasbMeanCon=tasbMeanCon[np.newaxis,:,:,:]

        if year==2019:
            PerturbedAll=taspMeanCon
            BaseAll=tasbMeanCon
        else:
            PerturbedAll=np.ma.concatenate((PerturbedAll,taspMeanCon),axis=0)
            BaseAll=np.ma.concatenate((BaseAll,tasbMeanCon),axis=0)
    # Five year mean
    PerturbedMean0=np.nanmean(PerturbedAll,axis=0)
    BaseMean0=np.nanmean(BaseAll,axis=0)
    # Here 96:192 means only use data in NH
    PerturbedMean=PerturbedMean0[:,96:192,:]
    BaseMean=BaseMean0[:,96:192,:]
    return PerturbedMean,BaseMean

def MainRead(Perturbedname,Basename):
    """
    MainRead
    Description: a main function to calculate duirnal cycles of mean tas, huv, and Bowen ratio 
                in the North hemisphere for the perturbed and base land tile by calling ReadVar and ReadVarNoMean 

    Input： Perturbedname--the name of the perturned land tiles (u_Urban/l_Lake/g_Grass)
           Basename--the name of the base land tiles (r_Rural/r_Rural/t_Tree)
    Output: tasp--duirnal cycles of mean air temperature in the North hemisphere for the perturbed land tile
    tasb--duirnal cycles of mean air temperature in the North hemisphere for the base land tile
    huvp--duirnal cycles of mean vapor pressure in the North hemisphere for the perturbed land tile
    huvb--duirnal cycles of mean vapor pressure in the North hemisphere for the base land tile
    PerturbedBowenRatio--duirnal cycles of mean Bowen ratio in the North hemisphere for the perturbed land tile
    BaseBowenRatio--duirnal cycles of mean Bowen ratio in the North hemisphere for the base land tile
    """ 
    [tasp,tasb]=ReadVar('tas',Perturbedname,Basename)
    [huvp,huvb]=ReadVar('huv',Perturbedname,Basename)
    
    [hflsuNoMean,hflsrNoMean]=ReadVarNoMean('hfls',Perturbedname,Basename)
    [hfssuNoMean,hfssrNoMean]=ReadVarNoMean('hfss',Perturbedname,Basename)
    #BowenRatio
    BowenRatioPerturbed = hfssuNoMean/hflsuNoMean
    BowenRatioBase = hfssrNoMean/hflsrNoMean
    # When the sensible and latent heat flux are too small (<0.15 W/m2) or the Bowen ratio is too larger(>25), the result is not quite meaningful
    PerturbedBowenRatio = np.ma.masked_where((np.absolute(hflsuNoMean)<0.15)|(np.absolute(hfssuNoMean)<0.15)|(np.absolute(BowenRatioPerturbed)>25),BowenRatioPerturbed)
    BaseBowenRatio =np.ma.masked_where((np.absolute(hflsrNoMean)<0.15)|(np.absolute(hfssrNoMean)<0.15)|(np.absolute(BowenRatioBase)>25),BowenRatioBase) #
    # Make sure the unmasked grids are consistent for the perturbed and base land tile
    PerturbedBowenRatio = np.ma.masked_where(np.ma.getmask(BaseBowenRatio),PerturbedBowenRatio)
    BaseBowenRatio= np.ma.masked_where(np.ma.getmask(PerturbedBowenRatio),BaseBowenRatio)
    
    PerturbedBowenRatio=AreaWeightedMean(PerturbedBowenRatio)
    BaseBowenRatio=AreaWeightedMean(BaseBowenRatio)
    return tasp,tasb,huvp,huvb,PerturbedBowenRatio,BaseBowenRatio


#----------------------------Urbanization-------------------------------
UrbannameInput='u_Urban';RuralnameInput='r_Rural';
# Read two arrays from the perturbed and base land tiles
# Their maskes are useful to ensure only grids with both perturbed and base land tiles are used
MaskU=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+UrbannameInput+'_Hourly_2019.nc')
MaskR=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+RuralnameInput+'_Hourly_2019.nc')
MaskArray=MaskU.variables['tas'+UrbannameInput[0]][0:24,:,:]-MaskR.variables['tas'+RuralnameInput[0]][0:24,:,:]
[tasUrban,tasRural,huvUrban,huvRural,BRUrban,BRRural]=MainRead(UrbannameInput,RuralnameInput)
del MaskArray

#--------------------------Lake-----------------------------
UrbannameInput='l_Lake';RuralnameInput='r_Rural';
# Read two arrays from the perturbed and base land tiles
# Their maskes are useful to ensure only grids with both perturbed and base land tiles are used
MaskU=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+UrbannameInput+'_Hourly_2019.nc')
MaskR=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+RuralnameInput+'_Hourly_2019.nc')
MaskArray=MaskU.variables['tas'+UrbannameInput[0]][0:24,:,:]-MaskR.variables['tas'+RuralnameInput[0]][0:24,:,:]
[tasLake,tasNonlake,huvLake,huvNonlake,BRLake,BRNonlake]=MainRead(UrbannameInput,RuralnameInput)
del MaskArray

#--------------------------Forest-----------------------------
UrbannameInput='g_Grass';RuralnameInput='t_Tree'
# Read two arrays from the perturbed and base land tiles
# Their maskes are useful to ensure only grids with both perturbed and base land tiles are used
MaskU=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+UrbannameInput+'_Hourly_2019.nc')
MaskR=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+RuralnameInput+'_Hourly_2019.nc')
MaskArray=MaskU.variables['tas'+UrbannameInput[0]][0:24,:,:]-MaskR.variables['tas'+RuralnameInput[0]][0:24,:,:]
[tasGrass,tasTree,huvGrass,huvTree,BRGrass,BRTree]=MainRead(UrbannameInput,RuralnameInput)
del MaskArray

# Define colors for three land tile pairs
UrbanColor='#CC0000';RuralColor='#16A085'
GrassColor='sandybrown';TreeColor='forestgreen'
LakeColor='cornflowerblue';NonLakeColor='salmon'

fig,ax = plt.subplots(nrows=3, ncols=3, figsize=(10,4.7))
Hours=np.arange(1,25,1)

ax[0,0].plot(Hours,tasUrban[:]-273.15,color=UrbanColor,linewidth=1.4,label="Urban Tas")
ax[0,0].plot(Hours,tasRural[:]-273.15,color=RuralColor,linewidth=1.4,label="Rural Tas")
# ax[0,0].set_xlabel('Local time',fontsize=14,labelpad=1)
# plt.ylabel('$\it{Ta}$'+' (°C)',fontsize=11,labelpad=1)
# ax[0,0].set_ylabel('tas (°C)',fontsize=11,labelpad=1)
ax[0,0].set_xticks(np.arange(0, 28, 4))
ax[0,0].set_xticklabels([ ]*len(np.arange(0, 28, 4)),fontsize=12)
ax[0,0].tick_params(axis="both", labelsize=11,direction='in')
ax[0,0].set_ylim( (12, 24) )
ax[0,0].set_xlim((-0.5, 24.5))
ax[0,0].set_yticks(np.arange(14, 23, 4))
deltax=ax[0,0].get_xlim()[1]-ax[0,0].get_xlim()[0]
deltay=ax[0,0].get_ylim()[1]-ax[0,0].get_ylim()[0]
ax[0,0].text(ax[0,0].get_xlim()[0]+deltax*0.03, ax[0,0].get_ylim()[1]-deltay*0.04, '(a)', horizontalalignment='left',verticalalignment='top', size=12)
ax0 = ax[0,0].twinx()
ax0.plot(Hours, tasUrban[:]-tasRural[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Urban-Rural Ta")
ax0.set_yticks(np.arange(0.3, 1.3, 0.3))
ax0.set_ylim( (0,1.4) )
ax0.tick_params(axis="y",direction='in')
# ax0.set_ylabel('Δtas (°C)',fontsize=9,labelpad=1)

# ax[0,1].plot(Hours,GridMeanVAPOR_PRES[:],'k',linewidth=1.4,label="Grid vapore pressure")
ax[0,1].plot(Hours,huvUrban[:]*0.01,color=UrbanColor,linewidth=1.4,label="Urban huv")
ax[0,1].plot(Hours,huvRural[:]*0.01,color=RuralColor,linewidth=1.4,label="Rural huv")
# ax[0,1].set_xlabel('Local time',fontsize=14,labelpad=1)
# ax[0,1].set_ylabel('$\it{Ea}$'+' (hPa)',fontsize=11,labelpad=1)
# ax[0,1].set_ylabel('huv (hPa)',fontsize=11,labelpad=1)
ax[0,1].set_xticks(np.arange(0, 28, 4))
ax[0,1].set_xticklabels([ ]*len(np.arange(0, 28, 4)),fontsize=12)
ax[0,1].tick_params(axis="both", labelsize=11,direction='in')
ax[0,1].set_ylim( (12.5, 15) )
ax[0,1].set_xlim((-0.5, 24.5))
deltax=ax[0,1].get_xlim()[1]-ax[0,1].get_xlim()[0]
deltay=ax[0,1].get_ylim()[1]-ax[0,1].get_ylim()[0]
ax[0,1].text(ax[0,1].get_xlim()[0]+deltax*0.03, ax[0,1].get_ylim()[1]-deltay*0.04, '(b)', horizontalalignment='left',verticalalignment='top', size=12)
# ax[0,1].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=2,loc='lower left', prop={'size': 5.7}) #, ncol=3
ax1 = ax[0,1].twinx()
ax1.plot(Hours, (huvUrban[:]-huvRural[:])*0.01,color='dimgrey',linestyle='dashed',linewidth=0.8,label="Urban-Rural huv")
ax1.tick_params(axis="y",direction='in')
ax1.set_ylim( (-0.9, 0) )
ax1.set_yticks(np.arange(-0.7, 0, 0.3))

# ax[0,1].plot(Hours,GridBowenRatio[:],'k',linewidth=1.4,label="Grid Bowen Ratio")
ax[0,2].plot(Hours,BRUrban[:],color=UrbanColor,linewidth=1.4,label="Urban")
ax[0,2].plot(Hours,BRRural[:],color=RuralColor,linewidth=1.4,label="Rural")
# ax[0,2].set_xlabel('Local time',fontsize=14,labelpad=1)
# ax[0,2].set_ylabel('Bowen Ratio',fontsize=11,labelpad=1)
ax[0,2].set_xticks(np.arange(0, 28, 4))
ax[0,2].set_xticklabels([ ]*len(np.arange(0, 28, 4)),fontsize=12)
ax[0,2].tick_params(axis="both", labelsize=11,direction='in')
ax[0,2].set_ylim( (-6, 6) )
ax[0,2].set_xlim((-0.5, 24.5))
deltax=ax[0,2].get_xlim()[1]-ax[0,2].get_xlim()[0]
deltay=ax[0,2].get_ylim()[1]-ax[0,2].get_ylim()[0]
ax[0,2].text(ax[0,2].get_xlim()[0]+deltax*0.03, ax[0,2].get_ylim()[1]-deltay*0.04, '(c)', horizontalalignment='left',verticalalignment='top', size=12)
# ax[0,2].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=2,loc='lower left', prop={'size': 5.7}) #, ncol=3
ax2 = ax[0,2].twinx()
ax2.plot(Hours, BRUrban[:]-BRRural[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Urban - Rural")
ax2.tick_params(axis="y",direction='in')
ax2.set_ylim( (0.7, 3.6) )
ax2.set_yticks(np.arange(1, 4, 1))
lines, labels = ax[0,2].get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2,bbox_to_anchor=(1.15,0.24),frameon=False, ncol=1,loc='lower left', prop={'size': 7})

# ax2.set_ylabel('ΔBowen Ratio',fontsize=9,labelpad=1)

ax[1,0].plot(Hours,tasGrass[:]-273.15,color=GrassColor,linewidth=1.4,label="Grass Tas")
ax[1,0].plot(Hours,tasTree[:]-273.15,color=TreeColor,linewidth=1.4,label="Tree Tas")
# ax[1,0].set_xlabel('Local time',fontsize=14,labelpad=1)
# plt.ylabel('$\it{Ta}$'+' (°C)',fontsize=11,labelpad=1)
ax[1,0].set_ylabel('tas (°C)',fontsize=13,labelpad=3,fontweight='bold')
ax[1,0].set_xticks(np.arange(0, 28, 4))
ax[1,0].set_xticklabels([ ]*len(np.arange(0, 28, 4)),fontsize=12)
ax[1,0].tick_params(axis="both", labelsize=11,direction='in')
ax[1,0].set_xlim((-0.5, 24.5))
ax[1,0].set_ylim( (8, 17) )
ax[1,0].set_yticks(np.arange(10, 16, 2.5))
deltax=ax[1,0].get_xlim()[1]-ax[1,0].get_xlim()[0]
deltay=ax[1,0].get_ylim()[1]-ax[1,0].get_ylim()[0]
ax[1,0].text(ax[1,0].get_xlim()[0]+deltax*0.03, ax[1,0].get_ylim()[1]-deltay*0.04, '(d)', horizontalalignment='left',verticalalignment='top', size=12)
# ax[1,0].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=3,loc='lower left', prop={'size': 6.5}) #

ax0 = ax[1,0].twinx()
ax0.plot(Hours, tasGrass[:]-tasTree[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Grass-Tree Ta")
ax0.set_yticks(np.arange(-0.3, 0.65, 0.3))
ax0.set_ylim( (-0.45, 0.65) )
ax0.tick_params(axis="y",direction='in')

# ax0.set_ylabel('Δtas (°C)',fontsize=9,labelpad=1)

# ax[1,1].plot(Hours,GridMeanVAPOR_PRES[:],'k',linewidth=1.4,label="Grid vapore pressure")
ax[1,1].plot(Hours,huvGrass[:]*0.01,color=GrassColor,linewidth=1.4,label="Grass huv")
ax[1,1].plot(Hours,huvTree[:]*0.01,color=TreeColor,linewidth=1.4,label="Tree huv")
# ax[1,1].set_xlabel('Local time',fontsize=14,labelpad=1)
# ax[1,1].set_ylabel('$\it{Ea}$'+' (hPa)',fontsize=11,labelpad=1)
ax[1,1].set_ylabel('huv (hPa)',fontsize=13,labelpad=0,fontweight='bold')
ax[1,1].set_xticks(np.arange(0, 28, 4))
ax[1,1].set_xticklabels([ ]*len(np.arange(0, 28, 4)),fontsize=12)
ax[1,1].tick_params(axis="both", labelsize=11,direction='in')
ax[1,1].set_xlim((-0.5, 24.5))
ax[1,1].set_ylim( (10.5, 12) )
ax[1,1].set_yticks(np.arange(10.8, 12, 0.5))
deltax=ax[1,1].get_xlim()[1]-ax[1,1].get_xlim()[0]
deltay=ax[1,1].get_ylim()[1]-ax[1,1].get_ylim()[0]
ax[1,1].text(ax[1,1].get_xlim()[0]+deltax*0.03, ax[1,1].get_ylim()[1]-deltay*0.04, '(e)', horizontalalignment='left',verticalalignment='top', size=12)

ax1 = ax[1,1].twinx()
ax1.plot(Hours, (huvGrass[:]-huvTree[:])*0.01,color='dimgrey',linestyle='dashed',linewidth=0.8,label="Grass-Tree huv")
ax1.set_yticks(np.arange(-0.4, 0.25, 0.3))
ax1.tick_params(axis="y",direction='in')
ax1.set_yticks(np.arange(-0.3, 0.2, 0.2))
# ax1.set_ylabel('Δhuv (hPa)',fontsize=9,labelpad=1)

# ax[1,1].plot(Hours,GridBowenRatio[:],'k',linewidth=1.4,label="Grid Bowen Ratio")
ax[1,2].plot(Hours,BRGrass[:],color=GrassColor,linewidth=1.4,label="Grass")
ax[1,2].plot(Hours,BRTree[:],color=TreeColor,linewidth=1.4,label="Tree")
# ax[1,2].set_xlabel('Local time',fontsize=14,labelpad=1)
ax[1,2].set_ylabel('Bowen ratio',fontsize=12,labelpad=3,fontweight='bold')
ax[1,2].set_xticks(np.arange(0, 28, 4))
ax[1,2].set_xticklabels([ ]*len(np.arange(0, 28, 4)),fontsize=12)
ax[1,2].tick_params(axis="both", labelsize=11,direction='in')
ax[1,2].set_xlim((-0.5, 24.5))
ax[1,2].set_ylim( (-6, 6) )
deltax=ax[1,2].get_xlim()[1]-ax[1,2].get_xlim()[0]
deltay=ax[1,2].get_ylim()[1]-ax[1,2].get_ylim()[0]
ax[1,2].text(ax[1,2].get_xlim()[0]+deltax*0.03, ax[1,2].get_ylim()[1]-deltay*0.04, '(f)', horizontalalignment='left',verticalalignment='top', size=12)

# ax[1,2].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=2,loc='lower left', prop={'size': 5.7}) #, ncol=3
ax2 = ax[1,2].twinx()
ax2.plot(Hours, BRGrass[:]-BRTree[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Grass-Tree")
ax2.tick_params(axis="y",direction='in')
ax2.set_ylim( (-1.7, 1.5) )
ax2.set_yticks(np.arange(-1, 1.2, 1))

lines, labels = ax[1,2].get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2,bbox_to_anchor=(1.15,0.24),frameon=False, ncol=1,loc='lower left', prop={'size': 7})

ax[2,0].plot(Hours,tasLake[:]-273.15,color=LakeColor,linewidth=1.4,label="Lake Tas")
ax[2,0].plot(Hours,tasNonlake[:]-273.15,color=NonLakeColor,linewidth=1.4,label="Nonlake Tas")
ax[2,0].set_xlabel('Local time',fontsize=14,labelpad=1)
# plt.ylabel('$\it{Ta}$'+' (°C)',fontsize=11,labelpad=1)
# ax[2,0].set_ylabel('tas (°C)',fontsize=11,labelpad=1)
ax[2,0].set_xticks(np.arange(0, 28, 4))
# ax[2,0].set_xticklabels([ ]*len(np.arange(0, 28, 4)),fontsize=12)
ax[2,0].tick_params(axis="both", labelsize=11,direction='in')
ax[2,0].set_ylim((4.5,12))
ax[2,0].set_yticks(np.arange(5, 10.2, 2.5))
ax[2,0].set_xlim((-0.5, 24.5))
deltax=ax[2,0].get_xlim()[1]-ax[2,0].get_xlim()[0]
deltay=ax[2,0].get_ylim()[1]-ax[2,0].get_ylim()[0]
ax[2,0].text(ax[2,0].get_xlim()[0]+deltax*0.03, ax[2,0].get_ylim()[1]-deltay*0.04, '(g)', horizontalalignment='left',verticalalignment='top', size=12)

# ax[2,0].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=3,loc='lower left', prop={'size': 6.5}) #

ax0 = ax[2,0].twinx()
ax0.plot(Hours, tasLake[:]-tasNonlake[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Lake-Nonlake Ta")
ax0.tick_params(axis="y",direction='in')
ax0.set_ylim((-1.4,1.3))
ax0.set_yticks(np.arange(-1.2, 1, 1))

# ax[2,1].plot(Hours,GridMeanVAPOR_PRES[:],'k',linewidth=1.4,label="Grid vapore pressure")
ax[2,1].plot(Hours,huvLake[:]*0.01,color=LakeColor,linewidth=1.4,label="Lake huv")
ax[2,1].plot(Hours,huvNonlake[:]*0.01,color=NonLakeColor,linewidth=1.4,label="Nonlake huv")
ax[2,1].set_xlabel('Local time',fontsize=14,labelpad=1)
# ax[2,1].set_ylabel('$\it{Ea}$'+' (hPa)',fontsize=11,labelpad=1)
# ax[2,1].set_ylabel('huv (hPa)',fontsize=11,labelpad=1)
ax[2,1].set_xticks(np.arange(0, 28, 4))
# ax[2,1].set_xticklabels([ ]*len(np.arange(0, 28, 4)),fontsize=12)
ax[2,1].tick_params(axis="both", labelsize=11,direction='in')
ax[2,1].set_ylim( (9.7, 11.5) )
ax[2,1].set_xlim((-0.5, 24.5))
ax[2,1].set_yticks(np.arange(10.2, 11.5, 0.5))
deltax=ax[2,1].get_xlim()[1]-ax[2,1].get_xlim()[0]
deltay=ax[2,1].get_ylim()[1]-ax[2,1].get_ylim()[0]
ax[2,1].text(ax[2,1].get_xlim()[0]+deltax*0.03, ax[2,1].get_ylim()[1]-deltay*0.04, '(h)', horizontalalignment='left',verticalalignment='top', size=12)

# ax[2,1].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=2,loc='lower left', prop={'size': 5.7}) #, ncol=3
ax1 = ax[2,1].twinx()
ax1.plot(Hours, (huvLake[:]-huvNonlake[:])*0.01,color='dimgrey',linestyle='dashed',linewidth=0.8,label="Lake-Nonlake huv")
ax1.tick_params(axis="y",direction='in')
ax1.set_ylim( (0.2, 1.4) )
ax1.set_yticks(np.arange(0.4, 1.3, 0.4))

# ax1.set_ylabel('Δhuv (hPa)',fontsize=9,labelpad=1)

# ax[2,1].plot(Hours,GridBowenRatio[:],'k',linewidth=1.4,label="Grid Bowen Ratio")
ax[2,2].plot(Hours,BRLake[:],color=LakeColor,linewidth=1.4,label="Lake")
ax[2,2].plot(Hours,BRNonlake[:],color=NonLakeColor,linewidth=1.4,label="Nonlake")
ax[2,2].set_xlabel('Local time',fontsize=14,labelpad=1)
# ax[2,2].set_ylabel('Bowen Ratio',fontsize=11,labelpad=1)
ax[2,2].set_xticks(np.arange(0, 28, 4))
# ax[2,2].set_xticklabels([ ]*len(np.arange(0, 28, 4)),fontsize=12)
ax[2,2].tick_params(axis="both", labelsize=11,direction='in')
ax[2,2].set_ylim( (-6, 6) )
ax[2,2].set_xlim((-0.5, 24.5))
deltax=ax[2,2].get_xlim()[1]-ax[2,2].get_xlim()[0]
deltay=ax[2,2].get_ylim()[1]-ax[2,2].get_ylim()[0]
ax[2,2].text(ax[2,2].get_xlim()[0]+deltax*0.03, ax[2,2].get_ylim()[1]-deltay*0.04, '(i)', horizontalalignment='left',verticalalignment='top', size=12)

# ax[2,2].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=2,loc='lower left', prop={'size': 5.7}) #, ncol=3
ax2 = ax[2,2].twinx()
ax2.plot(Hours, BRLake[:]-BRNonlake[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Lake-Nonlake")
ax2.tick_params(axis="y",direction='in')
ax2.set_ylim( (-2.5, 6) )
lines, labels = ax[2,2].get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2,bbox_to_anchor=(1.15,0.24),frameon=False, ncol=1,loc='lower left', prop={'size': 7})
 
plt.subplots_adjust(top=0.95,
bottom=0.19,
left=0.096,
right=0.805,
hspace=0.055,
wspace=0.605) 

#Save Figure
# plt.savefig('E:\\Keer_work\\DataPaper\\Figure5.png', dpi=450)
# plt.close()

